
%calibration
%estimating alp
cali=@(x)sum(grpsize.*( (ex_ante_allocation(x,updatek,ordering_list)-group).^2 ),'all'); 
x0=2*ones(size(group(1:end)));
A = [];
b = [];
Aeq = [];
beq = [];
nonlcon=[];
lb= zeros(size(x0));
ub=[];
options = optimoptions(@fmincon,'Display','iter','Algorithm','sqp','MaxFunctionEvaluations',1300,'MaxIterations',100,'OptimalityTolerance',1e-5,'StepTolerance',1e-4,'UseParallel',true);
%
tic
[alp,distance_po]= fmincon(cali,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
toc
% 
alp_po=alp;%[alp;5];
grpalp=alp;%[alp;5];
% % % % 
% % % calculating r^*
worstofffun=@(type)interim_total(alp_po,type,ordering_list,grpsize)-1;
worstoff=fzero(worstofffun,[0.001,1]);
r_ast_po=zeros(1,ngroup);
parfor agent=1:ngroup
     comb=ordering_list{agent}{1};
     lower=ordering_list{agent}{2};
     q_agent=ordering_list{agent}{3};
     numberpossibilities=ordering_list{agent}{4};
     r_ast_po(agent)=interim_allocation(q_agent,grpalp,worstoff,comb,lower,numberpossibilities); 
end
r_ast_vec_po=r_ast_po(ikc); %1*570 
r_fac_po=r_ast_vec_po(ia2);
alp_fac_po=alp_po(ikc);
alp_fac_po=alp_fac_po(ia2);
r_ast_po;


function[eq]=ex_ante_allocation(grpalp,grpk,ordering_list)
    nagent=length(grpk);
    eq=zeros(nagent,1);
    for agent=1:nagent
        a=grpalp(agent);
        comb=ordering_list{agent}{1};
        lower=ordering_list{agent}{2};
        q_agent=ordering_list{agent}{3};
        numberpossibilities=ordering_list{agent}{4};
        fun=@(type)density(a,type).*interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        eq(agent)=integral(fun,0,1);
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
    end
end
function[eq]=info_rent(grpalp,grpk,ordering_list)
    nagent=length(grpk);
    eq=zeros(nagent,1);
    for agent=1:nagent
        a=grpalp(agent);
        comb=ordering_list{agent}{1};
        lower=ordering_list{agent}{2};
        q_agent=ordering_list{agent}{3};
        numberpossibilities=ordering_list{agent}{4};
        fun=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        eq(agent)=integral(fun,0.01,0.99);
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
    end
end

function[suma]=interim_total(grpalp,type,ordering_list,grpsize) %remember to change interim total 
    nagent=length(grpalp);
    suma=0;
    parfor j=1:nagent
        comb=ordering_list{j}{1};
        lower=ordering_list{j}{2};
        q_agent=ordering_list{j}{3};
        numberpossibilities=ordering_list{j}{4};
        suma=suma+grpsize(j)*interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
    end    
end

%v1:
% function[agent]=interim_allocation(grp,grpalp,type,grpk,comb,lower) %alp needs to be a col vector   
%     q_agent=agent_allocation(grp,comb,grpk);
%     alp=grpalp;
%     higher=comb;
%     ntype=length(type);
%     agent=zeros(size(type));
%     parfor tt=1:ntype
%         prob11=prob1(alp,type(tt),higher);
%         prob22=prob2(alp,type(tt),lower);
%         agent(tt)=sum(prob11.*prob22.*q_agent,'all');
%     end
% end

%v2: use v2 of ordering_list
function[agent]=interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities) %alp needs to be a col vector   
    %q_agent=agent_allocation(grp,comb,grpk);
    alp=grpalp;
    higher=comb;
    ntype=length(type);
    agent=zeros(size(type));
    parfor tt=1:ntype
        prob11=prob1(alp,type(tt),higher);
        prob22=prob2(alp,type(tt),lower);
        agent(tt)=sum(prob11.*prob22.*q_agent.*numberpossibilities,'all');
    end
end

function[q]=realized_allocation(order,fac_updatek)
    nagent=length(order);
    q=zeros(nagent,1);
    rest=1;
    allocated=0;
    k=fac_updatek;
    for rk=1:nagent
        agent=order(rk);
        q(agent)=min(rest,k(agent));
        allocated=allocated+q(agent);
        rest=max(0,1-allocated);
    end
end


function[z1]=density(a,type)
    z1= a.*(1-type).^(a-1);
end

% function[z1]=buyertypetimesdensity(a,type)
%     z1= a.*(type).*(1-type).^(a-1)-(1-type).^a;
% end
% 
% function[z1]=sellertypetimesdensity(a,type)
%     z1= a.*(type).*(1-type).^(a-1)+(1-(1-type).^a);
% end

% function[t]=invdist(u,a) %given quantile, output type
%     t=  1-(1-u)^(1/a);
% end
function[z2]=prob1(a,type,higher)
    z2=prod(( (1-type).^a ).^higher,1)';
end

function[z3]=prob2(a,type,lower)
    z3=prod((1-(1-type).^a).^lower,1)';
end


